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Abstract. We consider generalized linear transient convection-diffusion prob- 
lems for differential forms on bounded domains in 1". These involve Lie 
derivatives with respect to a prescribed smooth vector field. We construct 
both new Eulerian and semi-Lagrangian approaches to the discretization of 
the Lie derivatives in the context of a Galerkin approximation based on dis- 
crete differential forms. Details of implementation are discussed as well as an 
application to the discretization of eddy current equations in moving media. 

1. Introduction 

Recall the classical linear transient 2nd-order convection-diffusion problem for 
an unknown scalar function u = u(x, f) on a bounded domain f2 C R n : 

dtu — eAu + (3 ■ grad u = f in Q , 

(1) u = on dfl , 

u(0) = u a . 

Here, and in the remainder of the paper, /3 : i— ► M. n stands for a smooth vector 
field. For < e <C 1 we encounter a singularly perturbed boundary value problem, 
whose stable discretization has attracted immense attention in numerical analysis, 
see [3 2) and the many references cited therein. 

The boundary value problem ([lj turns out to be a member of a larger family of 
2nd-order boundary value problems (BVP), which can conveniently be described 
using the calculus of differential forms. It permits us to state the generalized 
transient 2nd-order convection-diffusion problems as 

*d t w(t) - e{-l) l d * dui{t) + *Lpu{t) — ip in f2 C ffi™ , 

(2) i*lu = on dn , 

w(0) = loq . 

These are BVPs for an unknown time dependent /-form cu(t), < I < n, on the 
domain fi. The symbol * stands for a Hodge operator mapping an Z-form to an 
n — /-form, and d denotes the exterior derivative. Dirichlct boundary conditions arc 
imposed via the trace i*w of the Z-form uj. We refer to [2UJ Sect. 2] and [U Sect. 2] 
for more details and an introduction to the calculus of differential forms in the 
spirit of this article. By Lp we denote the Lie derivative of ui for the given velocity 
field (3 [2 Ch. 5] . It provides the convective part of the differential operator in @ . 
Details will be given in Sect. l2,ll below. 

In two and three dimensions differential forms can be modelled by means of 
functions and vector fields through so-called vector proxies, see [71 Sect. 7], |H 
Table 2.1], and [HI Table 2.10. Thus, for n = 3, in the case of * induced by 



Date: January 7, 2010, submitted to special issue of Disc. Cont. Dyn. Sys. 

Authors' affiliation: SAM, ETH Zurich, CH-8092 Zurich, 

{hiptmair,hcumann}@sam.math.ethz.ch. 

^Occasionally we will use the operator v. p. to indicate that a form is mapped to its corre- 
sponding Euclidean vector proxy 

1 



2 



HOLGER HEUMANN AND RALF HIPTMAIR 



the Euclidean metric on M 3 , the following vector analytic incarnations of ([2]) arc 
obtained. The solution Z-form w is replaced with an unknown vector field u or 
function u, cf. [19] Sect. 2]. 

• For I — we recover (JXJ) . 

• In the case I = 1 we arrive at a convection-diffusion problem 

<9fU + e curl curl u — (3 x curl u + grad(u • (3) = f in fi , 

(3) u x n = on dQ , 

u(0) = u . 

• The corresponding boundary value problems for vector proxies of 2-forms 
read 

dtu — e grad divu + /3(divu) + curl(u x (3) = f in , 

(4) u • n = on dQ , 

u(0) = u . 

• For / = 3 the diffusion term vanishes and we end up with pure convection 
problems for a scalar density u 

d t u + div(/3u) = / infi, 
() u(-,0) = u , 

where we assume (3 ■ n = on dVt, 

Though equivalent to ^ the vector proxy formulations very much conceal the 
common structure inherent in ©-((I]). Thus, in this article, we consistently adopt 
the perspective of differential forms. Thinking in terms of co-ordinate free differ- 
ential forms offers considerable benefits as regards the construction of structure 
preserving spatial discretizations. By now this is widely appreciated for boundary 
value problems for d * duj, where the so-called discrete exterior calculus [2"l [T3"ll2"2~] . 
or, equivalently, the mimetic finite difference approach 9,23 25 , or discrete Hodge- 
operators [8, 19] have shed new light on existing discretizations and paved the way 
for new numerical methods. 

All these methods have in common that Z-forms are approximated by Z-co- 
chains on generalized triangulations of the computational domain Q. This pre- 
serves deRham co-homology and, thus, plenty of algebraic properties enjoyed by 
the solutions of the continuous BVP carry over to the discrete setting, see, e.g., [SQl 
Sect. 2.2]. Moreover, simplicial and tensor product triangulations allow the exten- 
sion of co-chains to discrete differential forms, which furnishes structure preserving 
finite element methods. Also in this paper the focus will be on Galerkin discretiza- 
tion by means of discrete differential forms and how it can be used to deal with Lie 
derivatives. 

In light of the success of discrete differential forms, it seems worthwhile exploring 
their use for the more general equations © . This was pioneered in 6 , 23] and 
this article continues and extends these considerations. Revealing structure is not 
our only motivation: as already pointed out the discretization of ([1]) has been the 
object of intense research, but Q and (HJ) have received scant attention, though, ((3]) 
is relevant for numerical modelling, e.g., in magnetohydrodynamics. Apart from 
standard Galerkin finite element methods and [1, 29 , we would like to mention 
[TUGin] as attempts to devise a meaningful discretization of the transport term 
in Q. The first paper relies on a Lagrangian approach along with a divergence 
preserving projection, which can be interpreted as an interpolation onto co-chains. 
The second article employs edge degrees of freedom (1-co-chains) combined with a 
special transport scheme. Both refrain from utilizing exterior calculus. 

We stress that this article is confined to the derivation and formulation of fully 
discrete schemes for @ on a single fixed triangulation of f2. Principal attention is 
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paid to the treatment of the convective terms. The discussion centers on structural 
and algorithmic aspects, whereas rigorous numerical analysis of stability and con- 
vergence analysis will be addressed in some forthcoming work, because it relies on 
certain technical results that must be established for each degree / of the form u> 
separately. 

The plan of the paper is as follows: We first review the definition of Lie deriva- 
tives and material derivatives for differential forms and give a short summary on 
discrete forms. Next, we present a semi-Lagrangian approximation procedure for 
material derivatives of arbitrary discrete /-forms. Afterwards we present an Eulerian 
treatment of convection of discrete forms. These two sections focus on the formu- 
lation of fully discrete schemes. They are subsequently elaborated for the eddy 
current model in moving media as an example for non-scalar convection-diffusion, 
see also |17j . Finally we present numerical experiments for convection of 1-forms in 
2 dimensions. They give empiric hints on convergence and stability. 

2. Preliminaries 

2.1. Lie derivatives and material derivatives of forms. We introduce a space- 
time domain Qt := O x [to,ii] where fl C R n is a bounded curvilinear polyhedron 
and [toj^i] C R. Given a vector field (3 : Qt R n and initial values (x,s) € Qt 
the associated flow X(x,s,t) : Cl x [to,ti] 2 >-> R" is defined through 

^-X(x, s, t) = (3(X(x, s, t), t), t e [t , h] 
at 

X(x, s, s) = x . 

For fixed (x,s) the solution X(x, s, •) is called the characteristic curve through 
(x,s). A unique solution of problem (jB]) exists, whenever f3 is continuous in Qt 
and Lipschitz continuous in Q for fixed t S [fo,ti]- In this case 

(6) X(X{x,s,t),t,s) = x, VieO 
hence X(-, t, s) is the inverse of X(-, s, t): 

(7) X(;t,s)oX(-,s,t) = id. 

For convenience we abbreviate X s , t (-) — X(-, s,t). It follows directly from ([5]) that 
the Jacobian DX Sit (x) solves: 

(8) j t DX s , t {x) = D(3(X S}t (x), t)DX s>t {x) . 

We write for the space of differential /-forms on f2 (cf. Chapter V, 

Section 3]). Differential /-forms can be seen as continuous additive mappings from 
the set <S'(f2) of compact oriented, piecewise smooth, /-dimensional sub- manifolds 
of fi into the real number^. The spaces can be equipped with the i 2 -norm 

IMIo : = In *^Aaj, see [H Sect. 2.2]. For Euclidean Hodge operator this agrees with 
the usual L 2 -norm of the vector proxies. Although it is very common to regard the 
mapping w : S l (il) h- > C, u> £ as integration we adopt the notation < ui, Mi > 

of a duality pairing instead of J M ui, Mi £ S l (£l). Only for / = n we will use the 
integral symbol. 

0-forms can be identified with real valued functions. Recall the definition of the 
directional derivative for a smooth function / : f2 1— > R: 

OS.grad/XM) := H ffi^O^M . 

T->t T — t 



This notion of differential forms can be rigorously defined with tools from geometric measure 
theory 1 1 6 II 3 1 1 ■ Here we forgo these technical aspects and appeal to intuition, cf. | 20l Sect. 2.1] 
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The Lie derivative Lp of /-forms lu £ J- l (tt) generalizes this. Note that the point 
evaluation of 0-forms is replaced with evaluation on /-dimensional oriented sub- 
manifolds Mi £ S l (f2) . Thus the Lie derivative of a /-form ui is 

(9) < Law(t), Mi >:= hm . 

T— >t T — t 

In terms of the pullback X£ T defined by 

(10) < X* >t lj, Mi >:=< u, X t ,r{Mi) > 
we can also write 

X t * w — uj 

(11) L{jw{t) := lim ^ - . 

Following [B] the extrusion Extt lT (/3,Mi) = {X t . s (x) : t < s < r, x £ Mi} is the 
union of flux lines emerging at Mi running from t to r (Figure [T]) . The orientation 
of Mi induces an orientation of the extrusion Extt lT (/3, Mi) such that the boundary 
is: 

(12) dExt t , T ((3, Mi) = X ttT (Mi) - M t - Ext t , T ((3, dM t ) . 

Plugging this into the definition of the Lie derivative (fTTj) we get by means of Stokes' 
theorem 

T M , f ,. < uj,dExt UT (f3,Mi) > + < u,Extt, r (P,dMi) > 
< Latj(t),Mi > = hm '■ 

(13) r_< 

_ Hm < &u,Ext t , T {l3,Mi) > + < oj,Ext t , T (l3,dM{) > 

T->i T — t 

The contraction operator is defined as the limit of the dual of the extrusion: 

< u,Extt, T (P,Mi) > 



(14) < ipcj(t),Mi >:= lim 



t 



and we recover from (fl~3|) Cartan's formula [23 page 142, prop. 5.3] for the Lie 
derivative: 

(15) Lpu = ipdu + dipuo . 

For 0-forms the second term vanishes, for n-forms the first one. 
For time dependent /-forms u>(t) £ F 1 ^) the limit value: 

x;Mt) - w(i) 



(16) Dpujit) := lim 
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corresponding to ([lip is the rate of change of the action of Z-forms in moving media, 
hence a material derivative [J51 page 62]. We deduce: 



(17) 



Dgoj(t) = hm — ! ; h lim — ■ 

t->* r — i T->t t — t 

d 



Remark 2.1. Since the exterior derivative and the Lie derivative commute we 
have: 

(18) D p d = dD p . 

As a consequence closed forms remain closed when they are advected by the material 
derivative. If Dpuj{t) = and duj(0) — then this property is preserved: duj{t) = 
0,Vt. 

Remark 2.2. For 1-forms to G ^(il) with il C R 3 with vector proxy u = v.p.(w) 
formulas (| 1 T[) cm<i (|15|) give a general convective term 

v. Y>.(Lpuj) = —(3 x curl u + V(/3 • u) , 

while for 1-forms CO G J-" 2 (f2) with vector proxy u = v. p.(cj) £/ie convective term 
reads 

v. p^L^w) = /3(divu) + curl(u x /3) . 

For closed 1-forms lu € J 7 '^), f2 C R 3 with vector proxy u = v. p. uj, which satisfies 
divu = 0, we find 

v. p^Z^w) = d t u + curl(u x (3) . 

1.1. Whitney forms. We rely on a version of discrete exterior calculus spawned 
by using discrete differential forms. We equip fl with a simplicial triangulation 0^. 
The sets of subsimplices of dimension I are denoted by Si (f2ft) and have cardinality 
Ni := \S l h \. In words, S®(£lh) is the set of vertices of ilh, S^f^) the set of edges, 
etc. We write yV l {Vth) C J 71 ^) for the space of Whitney /-forms defined on fi/j, 
see Sect. 23] and [H Sect. 3]. From the finite element point of view a Whitney 
Z-form u>h G W l (£lh) is a linear combination 

v, 

(19) w h =526{(x)w i 

i=l 

of certain basis functions associated with subsimplices of ft /,.. The basis functions 
b\ have compact support and vanish on all cells that have trivial intersection with 
subsimplex s\ € S l h (£lh). Locally, e.g. in each cell s£ with s\ H s*j = s-, the basis 
functions b\ can be expressed in terms of the barycentric coordinate functions Xj 
and whose gradients d\. The barycentric coordinate functions Aj are associated 
with vertices s° of cell s£, e.g. s-flsj! = sj. If I = (Io, • • • ij) denotes the index set 
of vertices Sj ; ■ • ■ s ?; °f subsimplex s' then 

i i 

(20) &ii = E^ 1 )'^ A dA ^- 

The construction of these spaces ensures that the traces i*uJh on element boundaries 
are continuous. The evaluation maps (?-)i=i : F l ^ ^ N ', ^(w) : =< s' >, s- G S l h 
are the corresponding degrees of freedoms, which give rise to the usual nodal inter- 
polation operators II' : J~ l (£l) n- W l (fth), also called DeRham maps [20 , Sect. 3.1], 

N, 

(21) n^ = E^(* 

i=l 
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Remark 2.3. Discrete differential forms are available for more general triangula- 
tions |7| Sect. 25]. The Whitney forms feature piecewise linear vector proxies, but 
can be extended to schemes with higher polynomial degree \2(K Sect. 3.4]. For the 
sake of simplicity we are not going to discuss those. 



3. Direct and adjoint semi-Lagrangian Galerkin formulation 

There are two basically different strategies to discretize material derivatives. We 
either approximate the difference quotient (1161) directly or use the formulation (|17[) 
in terms of partial time derivative of forms and Lie derivative. The former policy 
is referred to as Lagrangian and in our setting of a single fixed mesh it is known as 
semi-Lagrangian. It is explored in this section. The latter is called Eulerian and 
reduces to a method of lines approach based on a spatial discretization of the Lie 
derivatives. Its investigation is postponed to Sect. H) 

We begin the study of the semi-Lagrangian method with a look at the pure 
transport problem: given some <p(t) £ F l (Sl) and oj(0) £ W l (Slh) find oj(t) £ •F'(n), 
such that 

(22) Dpu = ip . 

Here and below we assume in the following that on the boundary of SI the normal 
components of the vector field (3 vanish, hence •, •) = SI. This is a restrictive 

but very common assumption for semi-Lagrangian methods. Only a few semi- 
Lagrangian methods like ELLAM [TT] can handle a non-vanishing velocity field on 
the boundary. 

There are two different variational formulations for the transport problem (|22p . 
The direct variational formulation reads: given <p(t) £ J- 1 (SI) and initial data w(0) £ 
T l (Sl) find uj(t) £ P(St), such that 

r xt t uj(t) - w(t) r 

(23) / lim l ' T V ' — An = / p i Aij, Vn £ J- n ~ (SI) . 

Jn T-t J n 

The adjoint variational formulation for (|2"2"|) reads: given (p(t) £ J 71 (SI) and initial 
data uj(0) £ T l (Sl) find u(t) £ T l (Sl), such that 

f u(t) AX* f ri -u(t) An f , 

(24) / lim -±1 ^ ^— '- — / ip(t) A rj , VV7 £ J^~ l (SI) . 

./o r-yt At ./o 



Lemma 3.1. The direct \23\) and the adjoint \24\) variational formulations are 
equivalent. 

Proof. Let lo £ ^(Sl) and r\ £ J rn ~ 1 (Sl). Since X t}S (Sl) = S} and X Tyt o X t , T = id: 



X* tT u;(T)Av= / X* r cj(r)Ar7= / u) A X* t n , 

Jx Ttt (n) Jn 

and the equivalence follows. □ 
Corollary 3.2. Under the assumption X tjT (f2) = Sl,Vi, r: 

/ Lpuj Ar] = - I oj A Lpi] , w £ F l (Sl),n £ .F" _i (fi). 
Jn Jn 

Replacing the limits in (|2"3"|) and (j2"4")l with a finite difference quotient yields semi- 
discrete timestepping schemes: Given ip(t) £ J- l (Sl) and u>(t — At) £ J- l (Sl) find 
uj(t) £ P(Sl), such that 

(25) f ^ - X *^ Mt ~ M) A n = / A V , Vr, £ j r "~ l (S}) 
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and given ip(t) E P(Cl) and w(t - At) E find oj(t) E P(0,), such that 

f uj(t) Art — uj(t — At) A X* At f n f . 

(26) / t^±±i= / y>(t) A r) , y V EJ- n - l (Q). 

Restricting the semi-discrete formulation to the discrete spaces W l (Qh), we end up 
with the following direct and adjoint schemes: Given ip(t) E J~ l (Q) and u>(t — At) E 
W l (il h ) find uj h (t) E W l {fl h ), such that 

r uj h (t) - X* Af uj h (t~At) r 

(27) / * — - \f " A Vh = / Mt) A fj h , VrjH E W l (n h ) 



Jn At Jn 

and find 0Jh(t) E W l {£lh), such that 

*u) h (t) A rj h - *uj h (t - At) A X*_ Att f) h 



(28) 



At 

*ip(t)Afj h , Vv h EW l (n h ). 



= / *<4 

In 



A significant advantage of the semi-Lagrangian approach is the preservation of 
closeness (see Remark I2.ip . This property is important in many physical applica- 
tions. The adjoint semi-Lagrangian timestepping fulfils this in a weak sense. 

Definition 3.3. A form uj E J 71 ^) is weakly closed if f n *ui A dip — 0,\A/> E 
P '{{}). 

For if = the adjoint semi-Lagrangian timestepping boils down to 

*u; h (t)Ad^ h = / *uj h (t-At)AX;_ Att dtp h ViphEW 1 - 1 ^). 
', Jn 

As the exterior derivative and pullback commute we conclude 
*u)h(t) A di>h = / *u h (t — At) A diph , 



with iph = X*_ At t^h- Hence uJh(t) is weakly closed if uJh(t — At) is weakly closed. 

Remark 3.4. Another advantage of semi-Lagrangian methods is the straightfor- 
ward stability analysis for the homogeneous transport problem Dpu^ = 0. The 
corresponding direct scheme reads: given lu^t) E yV l (fth), find iOh{t) € W (fift) 

*u) h (t)Afj h = / *Xl t _ At u h (t- At) A fjh , W] h EW l {n h ). 
Jn 

Testing with uih we immediately get the stability estimate: 

K(*)||o < ||X t V A ^(i-Ai)|| 
< C\\u> h (t-At)\\ 

with C = C(f3,At) > 0. For il C K 3 and Euclidean Hodge the concrete constants 
can be deduced from vector proxy representations of the pullback \20l p. 245]. The 
stability estimates then read: 

||w/,(t)||o < sup \ det(DX t - Mit (jc))\*\\u h (t-At)\\ , u h €W°(n h ); 



IK(t)||o < sup (\(det(DX t _ M>t (x))\3 \\DX t ,t-At(x)\\) IK(t- At)|| , oj h E W\n h ); 
xen K ' 

||w h (t)||o < sup (\det(DX t ,t-At(x))\*\\DX t -At,t( x )\\) IK(t - Af)|| , oj h EW 2 (fl h ); 
xen K ' 

||w h (t)||o < sup \det(DX t ^ At (x))\^ \\u h (t- At)\\ , u h € W 3 (fi ft ), 
xen 
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where \\ ■ \\ is the Euclidean matrix norm. Stability will depend in general on the 
distortion effected by X t .t-At- For I = 1,2 even the usual assumption div/3 = 
does not guarantee stability. 

It is important to note that both the stability estimates and the preservation of 
weak closedness as previously stated assume exact evaluation of the bilinear forms 
J Q *X* t _ At w h (x, t - At) A f] h {x)dx and J n *uj h (x, t — At) A X*_ Att rj h (x)dx. Any 
non-trivial problem will require the use of additional approximations. We propose 
the following three approximation steps. 

(i) First we use some numerical integrator to determine approximations X t ^' (s°) 
to the exact evolutions X ijt /(s°) of vertices s°. Note that the direct scheme 
requires to solve © backward in time, while for the adjoint scheme we 
calculate forward in time. The simplest numerical integrator for this first 
step is the forward Euler scheme, which gives Xt,v (s°) = s\ + (t' — t)(3{s®). 

(ii) Then we approximate the evolution X tt t' (x) of arbitrary points x as the 
convex combination of the approximated evolutions of surrounding vertices: 



This yields a piecewise linear, continuous approximation X t ^'{x) of the 
exact evolution X tt t'(x). 

Note that for small time steps and Lipschitz continuous (3 the approx- 
imating flow X t> t' maps the simplicial triangulation fi/j onto a new sim- 
plicial triangulation fi^, such that each subsimplex s\ of flh has an image 
s\ =: X tl t'(sl) in Q h . 
(iii) Finally, we apply the nodal interpolation operators of discrete forms [201 
Sect. 3.3] to map the transported discrete forms X£ t _ At uih(t — At) and 
Xt_ At i rjh onto the space of discrete forms. 
Let us illustrate these three steps for the case of 0-forms and 1-forms: To deter- 
mine now the interpolant of a discrete transported 0-form X£ t ,uih,u)h G W (fife), by 
linearity it is enough to consider the basis functions. Since < TL°uj, >=< u>, s° > 
for all vertices s° and 0-forms w the matrix operator f , with entries 



maps the expansions coefficients of ujh to the coefficients of 11° X^ t ,uJh- This means 

that in each time step we not only need to determine the points Xt^'is®) but also 
the location within the triangulation O^. To find the element, in which X tjt '(s°) is 
located, we trace the path of the trajectory from one element to the next. Based 
on this data the matrix entries (f3T)f can be assembled element by element (see fig. 
[5] for the direct scheme). 

The advantage of the second approximation step, the linear interpolation, is 
elucidated by the treatment of 1-forms. The interpolation of a transported discrete 
1-form X^ t ,LUh is determined by the interpolation of transported basis forms and 
the condition < rFA^^u;/,, — X^/UJh, sj >— for all edges sj € S^. This defines 
a matrix ~P\ t ,, mapping the expansion coefficients of u>h to those of X£ t ,u)h- The 
matrix entries 



(29) 




(30) 



p% := l9(X ttt ,) = A,(AV( S °)) = A,(AV(s°)) 



(31) Pb :=lHXtt'bl)=<Kt>bl" s t>> 



are line integrals along the straight line from Xt >T (sj) to X tjT (s°), if s\ is the edge 
connecting the vertices s° and s° (see fig. |3]for the direct scheme). To determine 
the entries of the i-th row, we trace the straight line from X t ji(sj) to X^t'is^) and 
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Figure 2. To determine the location of X t ,t-At(si), we move 
backward along the trajectory X tr (s®) starting from s° and iden- 
tify the crossed elements sj, sS^ S3 and s°. In this case p® kl p®i and 
p\ m are the only non-zero entries in the i-th row of Ptt-At- 




Figure 3. The transported edge X t ,t-Ats} (black curved line) is 
approximated by a straight line X t: t-&ts\ (black dashed line). In 
the case depicted here all basis functions associated with edges s\, 
of elements s™, s£ and S3 1 yield a non-zero entry p\ v . 



calculate for each crossed element the line integrals for the attached basis functions. 
If e.g. the line crosses an element s n with edge s\, from point a to point b (see fig. 
2] for the direct scheme), then the element contribution to p\ it is: 



b\ = / b\, =\ J ,(a)\ k ,(b)-\ k ,(a)\ J ,(b), 

X tit ,(sjna») J[a,b] 

where s®, and si, are the terminal vertices of edge s\, . 



10 



HOLGER HEUMANN AND RALF HIPTMAIR 





Figure 4. The line [a, b) is the intersection of the approximation 
of the transported edge with element s n . 

Observe that we do not need to calculate additional approximate solutions of 
([6]) to determine the matrix entries (j31j) . Piecewise linear interpolation of the 
approximate evolution of vertices automatically gives approximate evolutions of all 
subsimplices. 

Although the approximate flow Xt,t' is n °t smooth, we can show that the fully 
discrete adjoint scheme preserve the weak closedness property. Is is easily estab- 
lished that exterior derivative and interpolation of approximate pullbacks com- 
mute [201 Sect. 3.2]. 

Lemma 3.5. Let IF : (ft) i— > W l (£lh) be the standard nodal interpolation oper- 
ator and Xf tT (x) := X t _ T (xi)\i(x), where the coefficients Xf >T are approximate 
solutions of the characteristic ODE ([6]), then we have: 



Finally we state the fully discrete Semi-Lagrangian timestepping schemes for the 
general convection-diffusion problem @: 

Definition 3.6. Given a temporal mesh = <o < ti < < • ■ • > the discrete semi- 
Lagrangian timestepping scheme for the direct variational formulation of convection- 
diffusion of differential forms is: Given uih(to) <E W l (fth) and cp € J- l (ft), find 
^h(tk),k > such that: 



The discrete Semi-Lagrangian timestepping scheme for the adjoint variational for- 
mulation of convection- diffusion of differential forms is: Given oJh(to) £ W l (flh) 
and (p £ J- l (fl), find uih(tk),k > such that: 



dtt l xi T u = Tv l+1 x; T duj , e w l (n h ) ■ 



(32) 




(33) 
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Remark 3.7. The proposed sequence of approximation steps allows for a compre- 
hensive stability and convergence analysis of fully discrete schemes. Stability of the 
fully discrete schemes in f2 = K 3 for example follows from remark \3~4\ We only need 
to control perturbations due to some numerical integrator, the linear interpolation 
and interpolation of forms. We will address this in some future work. 

Remark 3.8. The standard finite element approach for approximating bilinear 
forms like J n *X^ t ,uJh A rjh would be quadrature. But since the product of X[[ t ,uih 
and rjh is smooth only on intersections s"nl(^(s"), this entails finding all intersec- 
tions of an element s™ of the triangulation with any other element X t /. t (s") of the 
transported triangulation, and apply some quadrature there. We conclude that such 
an approach is even more complex than methods based in interpolation. Further, it 
is totally unclear how such an approximation affects stability and conservation of 
closedness. 

Remark 3.9. // we use in the scalar case I = a vertex based quadrature rule to 
approximate the element integrals J s „ n°X t * r a;^(r) A n, the resulting scheme corre- 
sponds to a Lagrange- Galerkin scheme from J28f using vertex based quadrature. 

Remark 3.10. Notice that in both (|32[) and (1331) the unknown ujh(tk+i) is obtained 
by solving a standard discrete diffusion problem arising from a symmetric positive 
definite bilinear form, which are amenable to fast iterative solvers [3,18,21]. In the 
words of [33], the semi-Lagrangian method is "solver-friendly" . 

4. EULERIAN APPROACH: DISCRETE LlE-DERIVATIVES 

The Eulerian approach treats the discretization in space and the discretization 
in time separately. We therefore consider the stationary version of the convection- 
diffusion equation @: 

, . —e(—l) l d * duj + *Lpuj = tp in S] C 1™ , 

^ ' i*uj = ondn. 

Since the approximation of the diffusion term is well-known |19) . we will discuss 
only the approximation of the convective terms in (|3"4")l . We will pursue two different 
ways to obtain spatial discretization for the Lie-derivative terms. 

4.1. Standard scheme. The standard Galerkin technique to approximate 
J n *LpUh A rjh, 0Jh> f]h € VV (0^) relies on elementwise quadrature. While evalua- 
tion of the part *ipduih A rjh is standard in the finite element context, we en- 
counter some difficulties for the term J n *dipujh A rjh, which occurs for I > 0. The 
term dip implicitly requires spatial derivatives that fail to be square integrable for 
discrete differential forms. For instance, for 1-forms uj £ .F'(f2) with smooth vector 
proxy u := v. p.(w) we have: 

(35) grad(/3 • u) = (Df3) T u + {DuffS 

(36) = (Df3fu - (Du) asym (3 + {Du) sym [3 , 

while for discrete 1-forms Wh £ W l (fth) with vector proxy on each element 
the symmetric part (fu/,) ssm := ^(Duh + (flu| 1 ) T ) vanishes. A scheme based on 
elementwise quadrature will always miss (Du) sym f3. Stated differently, the exterior 
derivative dipUh of the contraction of a discrete form is not defined in a strong sense. 
But it is defined in a weak sense: For any smooth n — I form 7 £ F n ~ l (Q) , 1* ~f = 0, 
and discrete differential form uih £ W l (flh) elementwise integration by parts yields: 

N n „ N n 



(37) (-1)' f ipco h Ad-f = J2 I d h ipio h A"/-J2 f 



ipu h A 7 , 
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where dh is the restriction of d to each element s". Since each face inside Q occurs 
twice in the right sum we can rewrite this term as sum of face integrals over certain 
jump terms [ipUJh] 1 ^ 1 ■ 



(38) (-1)' / ipw h /\dj = Y^ / dhipuh A 7 - X! / i^h] 1 

Jn i=i Js ? j=i Js " 1 

One is now tempted to use the right hand side of (1551) to define approximations for 
the convection terms. The difficulty here are the face integrals, which are not well- 
defined for 7 G W l (flh)- Motivated by finite volume approximations, we replace 
in the case 7 G W'(f2/i) in the face integrals 7 with a flux function 7 = 7(7 , 7 r ) 
depending on the values 7' and j r on adjacent elements. 

We use the characterization (fTTj) of the Lie derivative to determine such consis- 
tent flux functions. 

Definition 4.1. Let (3 be Lipschitz continuous, u),fj G J 7 (CI) and At > then 

(39) &Ar(w,»j) := / *Lp T Lu A fj 



j\i 

is the variational central finite difference, with 

L e u - 2Ar~ ■ 

We could have used one-sided finite differences as well. Next we show that 
the limit limA-r-s.0 ^ArK, fjh), u)h> Vh G W (fifr) exists also for discrete differential 
forms. 

Lemma 4.2. The limit b (uj h ,fjh) := lim AT ^o &Ar {^h, Vh), ^h, Vh G W l (Qh) exists, 
and we have the representation 

(40) 

n„ K-i l 

b (u) h) fj h ) = *Lp(j h Afj h - ^2 -"h) A [vt + Vh )) > 

i=i Js r j=i Js 7° 

with upwind and downwind traces uj^ and uj7 and the pointwise limit LpuJh(x) = 
liniAr-s-o Lp T uJh{x). J° C {l,...iV n _i} is the index set of faces, that are not on 
the boundary ofCt. It's cardinality is N°_ 1 . Here we used 

(41) f i 3 K A ^) : = A lim . ./ ^"^ Ai. 

Proof. Let us first note that for all a; g" 9s" the pointwise limit LpUh(x) :— 
limAr-i-o Lp T ojh (x) exists. 

We decompose the central difference into a sum of upwind and downwind finite 
differences: 

, , ~ x 1 / f X t,t+Ar UJ h -Wh f U>h- X tt-A T Uh _ 
bAr K, Vh) = g * ^ A % + ^ * — A Vh 

and look first at the limit for the upwind finite difference. We split the integration 
over the domain fl in a sum of integrals over patches Pij(Ar) :— s™ n -^tt_A T ( s i) 
with smooth integrand: 



/ 

Jn 



' : A fj h = V / *— — — A % := ^ iij(AT) 



and distinguish three cases, see Fig. [5] 
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• i £ j and s™ n s™ = 0: 

(42) lim L J At) = 

Ar^O 

• i £ j and sj'ns™ = s£ _1 : This implies first that s£ _1 (£ dtt. If elements s[' 
and share a face s^" 1 it is clear that Pjj(Ar) = Ext t -& T ,t{@, s 7 ^ 1 ) 4 
0(At 2 ) and we have 

,. f [ UJ h „ 

lim / * — — Ar/h = hm / * — — A % 

(43) ^ d ^(At) At ^oJ Extt _ A ^ s n-i } At 

Js l 

and 

,. f Xt t _ Ar uJh _ /" Xt t _£ T u>h 

hm / * — j A 77/j = lim / * — A 77/,, 

(44) ^°^(At) At Ar-fO^,^^,,--!) At 

= - / h(* u h Af ih) 

Js k 

by the definition of contraction. 

• i = j : Since by assumption limA T ^o *— — ?_h±—^i!?h j s bounded on Pj ) j(Ar) 
we get: 



,. f - X* t _ Ar uj h _ 
hm / * A rjh 

At^o /p. i All At 

(45) P " (Al) 

*£j8Cd/, A fjh ■ 

After summation over all elements we obtain: 

f u h - X* t _ Ar uj h _ 

hm / * Ar)h = 

Ar^oJ n At lh 

(46) Nn „ if. 



i=l ""i j=l " V 

for the upwind finite difference. A similar argument for the downwind finite differ- 
ence shows: 

hm / * A r) h = 

Ar^oJ n At 

(47) N n _ K-i . 

*Lpuj h f\f) h -Y^ / ^ ~ w ft ) A ■ 

i=l Js ? 3 = 1 Js "o 

j 

Combining these two results for upwind and downwind finite difference we get the 
assertion. □ 



Remark 4.3. The careful examination of the proof of Lemma \4-^\ shows that the 
existence of the limit can be shown for quite general approximation spaces for dif- 
ferential forms. We used only the continuity inside each element. 

Remark 4.4. For the sake of completeness we rewrite here the standard convection 
bilinear form for vector proxies u,v or u, v in the case C K 3 . n denotes the 
normal vector of a face. Discrete Whitney 0- forms are continuous. Hence in this 
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J X t _ A r,t(a) 



Figure 5. Illustration for proof of Lemma 1431 In the case of the 
upwind finite difference we have, that discrete /-forms Qh G W l , I > 
are discontinuous on edges of the blue triangulation ilh- Their 
pullbacks X£ t _ AT LUh are discontinuous on edges of the red dashed 
triangulation X t -AT.t(^h), the image of flh- 

case the face integrals in (|40[) vanish and we obtain the standard bilinear form for 
scalar convection: 



(48) 



bo{vh,rjh)~ / /3 • gradu(x)v(x)dx . 
Jn 



For discrete 1- forms we have both volume and face contributions: 



N 3 



(49) b (uj h , fj h ) ~ ^ / gradJ/3 • u) • v - (J3 x curl u) • vdx 

i=i 

/3-n+(u+-u-).(v++v-)^, 



i=i"V° 



and similar for discrete 2 forms: 



N 3 



(50) b (cj h ,f) h ) ~ ^ / curU(u x /3) • v + y9- vdivudx 

i=i " /s f 



(v+ + v-)eZS. 



for discrete 3-forms we get: 



N 3 



3 = 1 -I? 



(51) 6o(w*.%)~S ^h{uf3)vdx--J2 f3-n+{u+-u-){v++v-)dS 



Remark 4.5. For 3-forms in 3D the limit (|40|) corresponds to the non- stabilized 
discontinuous Galerkin method for the advection problem from f 101. In the lowest 
order case this is identical with the scheme for discrete 3-forms. The limit for the 
upwind finite difference yields the stabilized discontinuous Galerkin methods with 
upwind numerical flux. 

In practice we will use local quadrature rules to evaluate the volume and face 
integrals in 6o(w/ l , fjh)- This will not destroy the consistency if the quadrature rules 
are first order consistent. 
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In the definition of 6at from Def. I4.1l we may also appeal to Lemma IXT1 and shift 
the difference quotient onto the test function. This carries over to the limit and 
yields two different Eulerian timestepping methods for problem ([2]) that, in analogy 
to Sect.[3l we may call "direct" and "adjoint". 

Definition 4.6. The direct standard Eulerian timestepping scheme for convection- 
diffusion of differential forms is: given LUh(to) € W l (£lh) and ip g .F (fi), find 
Wh{tk), k > such that: 

frn\ f Uh{tk+l) — Uh{tk) . - , f , / , N . 

(52) / * — — - — - — - A rjh + / *duj h (t k+1 ) A dr/ 
Jn tk+i — ife Jn 

+ 60(^(^+1),^)= / *tpAf) h , vfj h ew l (n h )- 

Jn 

The adjoint standard Eulerian timestepping is: given Uh(to) G W (fi/,) and ip £ 
.F (Q), find tUh(tk), k > such that: 

(53) / * — — Ar/ h + / *du h (t k +i) A drjh 
Jn tk+i ~ tk Jn 

- b (fj h ,uj h (t k+1 )) = *p A fj, Vfj h e w l (n h ) ■ 
Jn 

4.2. Upwind scheme. Now we use the definition (fTTj) of Lie-derivatives as limit 
of a difference quotient to define a proper interpolation operator for Lie-derivatives 
of discrete forms. Let sj £ S l h be a subsimplex of the triangulation. Then both 
limit values, the limit from the upwind direction 

54) <L 3 U) h ,s\>:= lim — : , meWSl/, 

p At^o+ At 

and the limit from the downwind direction 

(55) <L+u h ,a >:= lim M+Ar * ^— = — , u>h e W 1 ^) 

p Ar-s-o+ At 

exist for Lipschitz continuous (3. The continuity of traces of discrete differential 
forms ensures that integrals are well-defined. The existence of upwind and down- 
wind limits then follows from the existence of the limits for locally linear f3, which 
in turn can be calculated explicitly by distinguishing different geometric situations. 
We skip the tedious details. 
Observe that in general 

(56) < L+u h ,s\ >^< Lpuj h , s\ > , 
e.g. for O-forms a short calculation shows 

(57) < LpU h , S 1 >=/3(s° i ) - grade,; ( s °) , 
and 

(58) < L+u h , S ° >= /3(s°) ■ grad w+( S °) , 

with upwind and downwind traces lo7 and Only for cases where u>h is contin- 
uous in a neighbourhood of subsimplex s\ upwind and downwind limit agree. 

Since (|54[) and (|55[) correspond to integral evaluation of the usual interpolation 
operators, we could either use the upwind or the downwind limit to define discrete 
Lie-derivatives. We will use the upwind limit for the direct formulation and the 
downwind limit for the adjoint formulation. This is motivated by taking the cue 
from semi-Lagrangian schemes. 
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Definition 4.7. Let (3 be Lipschitz continuous and ojh € W (Oft). Further let 
(s|)^[ be the subsimplices corresponding to the degrees of freedom of basis functions 
(fife then 

(59) Lp ih u h -=^< L pu, s\ > b\ 

i 

is the upwind interpolated discrete Lie-derivative. Analogously the downwind inter- 
polated Lie- derivative is 

(60) L+ ih u> h :=J2< L 0^4>b l i - 



It follows from the Cartan formula ([T5]) that 

(61) < LpU} h ,s\ >=< d(ipU h ) +i^(duj h ),s l i >, 
with 

(62) <i g u h .s\>:= lim < u h , Ext^t-AriP, s i) > ■ 

Remark 4.8. The exterior derivatives can be evaluated exactly for discrete differ- 
ential forms. Thus, the contraction operators are the only approximations. But for 
locally constant discrete forms these approximations are exact and one can prove 
consistency 

(63) | f *Lp Afj- [ *Lp h n l u> Ml l fj\ = 0(h), Wui, fj e J"'(0) , 

Jn Jn 

once one has interpolation estimates for discrete differential forms. This is outside 
the scope of this paper. 



Remark 4.9. We finally want to stress that it is important to treat the limits (|54|) 
and (|55p for I > as limits of co-chains rather than integrals with certain integrands 
that are defined as pointwise limits. First, these pointwise limits are not well- 
defined since discrete forms are in general discontinuous across element boundaries. 

Second, even for special cases where the pointwise limit limA T ->o ~ ^ t 7 T AT "' 1 (x) 

exist we can have 

(64) lim < ^- X tt-A^H A 

v ' Ar^0+ At 1 r 

< lim " h ~^'- ATt % i>, u*eW(n fc ). 

At->o+ At 

To understand this, we may consider the 2D example depicted in Fig. [3] with con- 
stant (3. We take ujh € W 1 such that < uih,s\ >— for all edges sj <= S 1 except the 
edge s\ between vertices and s§ and calculate the projection of the Lie-derivative 
onto edge s\ between vertices s° and Sg. Then it is clear that 

(65) lim — - X W-*T Uh ( x )=0, xesiu{s°} 
and therefor 

(66) < lim ^"^ ,4 >=0- 
v ' Ar^0+ At 2 
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On the other hand, one can easily show that: 
< uj h - Xl t _ AT w h ,sl > 

(67) = Al (s° 2 - Ar/3) A 2 (s° - ^/3) - Ai (s° - A 2 (s° - Ar/3) 

= Ar(-^= grad Ai • - grad Ai • (3) + 0(Ar 2 ) , 

Zience 

hm < ,s 2 >= -=— gradAi • (3 ^ . 

At^o+ At y/2 

Since an explicit calculation of the upwind limits (|54j) might be expensive for arbi- 
trary (3, we have to find consistent approximations. The preceding discussion shows 
that this must be done very carefully. Replacing the integration simply with some 
quadrature will not yield a consistent approximation. 



Figure 6. Example that we can not change the order of integra- 
tion and limit. 

We combine the upwind discretization with an implicit Euler method and state 
the implicit upwind methods for convection-diffusion for discrete differential forms. 

Definition 4.10. The implicit direct upwind Eulerian timestepping scheme for 
the variational formulation of convection-diffusion of differential forms is: Given 
Wh(to) & W'(fi/j) and if £ find Uh(tk), k > such that: 

(69) / * — — An + e / *duj h (t k+1 ) A dn 

Jn t k+i — £fc Jn 



+ I *L ph uj h {t k+l ) Afj= I *tpAfj, yfjew'{n 



in Jn 
The adjoint upwind Eulerian timestepping scheme for the adjoint variational for- 
mulation of convection- diffusion of differential forms is: Given Ldh(to) G W l (flh) 
and if € J- l (il), find O0h(tk),k > such that: 

(70) / * — K — ^ — Av + e / *du) h (t k+ i) A dr) 

Jn l k+i — Efe Jn 

*u) h (tk+i) A L± h fj = / *ifAij 1 VijeW l {n h ). 
in Jn 

In contrast to the semi-Lagrangian schemes (l32|) and (j33jl . in ((69|) and (|70|) we 

face a non-symmetric algebraic system in each time step, whose iterative solution 

can be challenging. We therefor consider a second kind of upwind Eulerian scheme, 

that treats the convection in an explicit fashion. 
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Definition 4.11. The direct semi-implicit upwind Eulerian timestepping scheme 
for the variational formulation of convection-diffusion of differential forms is: Given 
w/i(£o) G W (f2ft) and <p € .F (fi), find u>h(t k ), k > suc/i i/wrf: 

(71) / * — ^— ^ — An + e *du h (t k+1 ) A dry 

Jn tfc+i — tfc in 



J /3 

TTie adjoint upwind Eulerian timestepping scheme for the variational formulation 
of convection- diffusion of differential forms is: Given uih(to) <= VV arad y € 

(fi), /ind 0Jh{tk),k > sucft £/ia£: 

(72) / * — ^— ^ — A?? + £ / *du} h (t k+ i) A dr) 

Jn tk+i ~ tk Jq 



*Uh(tk) A Lt h rj = / *tpAf), Vfj€W l {n h ). 
Ja 

These semi-implicit schemes do not coincide with the usual method of lines 
approach. Nevertheless there are two heuristic arguments that could justify the 
semi-implicit schemes. Note first that for discrete 0-forms the Semi-Lagrangian 
and the semi-implicit Eulerian scheme are identical, if we use a time step length 
At = 0(h) and explicit Euler for solving the characteristic equation ([5]). The 
direct Semi-Lagrangian scheme (p?2"j) for discrete 0-forms was, written here for vector 
proxies u and v for ojh and rjh & W°(flh)- Find u such that 

u(t k+1 ,x) - U°u(t k ,X t tk (x)) 

while the semi-implicit Eulerian scheme was: Find u such that 

f u(tk+i,x) u(tk, x) v ^ x ^ x _|_ g /* g rat i u (£ fc+1)X ) . gradw(a;)dx 

is) ^fe+i ~ tk Jn 

+ / L /3h u(t k ,x)v(x)dx = / f(x)v(x)dx, VueW°(fi), 
Jn ' Jo 

with 

A 

Lp h u(t k ,x) = ^2f3{sf) • gradu(i fc ,s°)| Ti Ai(a;) . 

i=l 

On the other hand Xt k+ll t k (s°) = — (tfc+i — t k )/3(s®) and by the time step 
condition At = 0(h): 

n°u(t k ,X tk+utk (x)) = J2 «(**, 
i=i 

- (ife+i - t k )[3(s°) ■ gradu(i fc ,s°)| Ti A,;(x) , 

hence we have the identity: 

U°u(t k ,X tk+lttk (x)) = u(tjfe, a;) - (t fc+1 - t k )L p h u(t k ) . 

The second argument, that supports the definition of semi-implicit schemes applies 
to forms of any degree and follows from the identity: 

X t ,t>w = ui — (X t ,t>u - w) 

= uj-(t' -t)L p uj + 0((t' -i) 2 ), weJ 1 )!!). 
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Hence we can derive the semi-implicit schemes (j6"9")l and (|70l) from the Semi-Lagrangian 
schemes (|3"2")) and (153"]) in replacing e.g. 

U l X th+1> t h Uh(tk) ~ Wh(tfc) - (t k - t k+ i)Lp h u> h (t k ) ■ 

Remark 4.12. TTie analysis of a semi-implicit discontinuous Galerkin scheme for 
scalar non-linear convection- diffusion can be found in Ji^] /. For the linear case this 
scheme is modulo certain quadrature rules identical with the semi-implicit Euler 
scheme for 3-forms. 

Remark 4.13. There is also a very close link between both the semi-Lagrangian and 
semi-implicit Eulerian methods and arbitrary Lagrangian-Eulerian methods (ALE) 
{3 p. 322]. ALE methods are operator splitting methods, that treat the diffusion part 
in a Lagrangian and the convection part in an Eulerian fashion. After a certain 
number of Lagrangian iterations steps an Eulerian iteration step maps the mesh 
function on the distorted mesh back to the initial mesh. The semi-Lagrangian and 
semi-implicit Eulerian methods apply this mapping in each time step. In many 
problems these remap operations should preserve certain properties of the mesh 
function, e.g. volume or closedness. It is the discrete differential form approach, 
that naturally guaranties such properties for the methods presented here and for the 
ALE method in JS~U$ . 

5. Application: Magnetic convection-diffusion 

In this section we will derive two different convection-diffusion equations for the 
electromagnetic part of magnetohydrodynamics (MHD) models. The first one is an 
equation for the magnetic field h and requires the adjoint methods while the second 
one for the magnetic vector potential a can be solved via the direct methods. 

Commonly in MHD applications, one neglects the displacement current. This 
reduced model, called eddy current model, is a system of equations for the magnetic 
field h 6 J rl (fi), the electric field e e J rl (fi), the magnetic field density b € F 2 {Vt), 
the current density j G .F 2 (f2) and external source / £ J" 2 (f2): 



(73) de = -d t b in ft , 

(74) dh = j + f inn, 

(75) j = * a (e — ipb) in 51 , 

(76) * li h = b in fl. 



A uniformly positive conductivity a will be assumed in the sequel. 

To rewrite the eddy current model in terms of a material derivative we substitute 
e = e+ipb and add —ipdb to Faraday's law (l73l) . This leaves the solution unchanged, 
since db = 0. Hence we end up with the system: 



(77) de = -Dpb in fl , 

(78) dh = j + f in fi , 

(79) j = * CT g in il , 

(80) * fl h = b mil. 



Next we state the two different variational formulations, cf. [20j Sect. 2.3]. 

5.1. /i-based variational formulation. For simplicity we impose homogeneous 
electric boundary conditions on <9f2. Testing (fTT|) with h' € J rl (J7), integration by 
parts yields: 

(81) / e A dh' + [ D p b A h' = . 

Jn Jn 
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We eliminate e using Ohm's law and end up with the following variational formu- 
lation: Seek h € J" 1 ^), b € J 72 ^) such that: 

(82) f D p bhh'+ [ * a -x{dh- f)f\dti = Q, V/i'e? 1 ^). 
Jn Jn 

(83) f * M feA/i"- f bAh" = 0, V^e^fi). 
Jn Jn 

In order to eliminate b as well we use Corollary 13.21 and get 



(84) / Ah' — *„h A L^fe') + / (dh -f)Adh' = 0, W e ^(tt) 

or, equivalently 
(85) 

* u h(t) Ah' -* u h(t- At) AX? Ait h' f , 

lim i *~ AM + / *„.-! {dh-f)Adh' = , V/i'eJ 1 

At At Jn 

This variational formulation can be approximated using a Semi-Lagrangian frame- 
work (|3"3"]) yielding an algebraic system: 

(86) (M M + AtC,-! )h< = (P ( 1 _ Ai( ) T M Al h t - At + Ati , 

where M M and C CT -i are the mass matrices J n *^cJ/j A % and f n * a -idu)h A dfjh, 
Vh: Vh S W 1 (il/j). The evaluation of the right hand side requires the calculation of 
the matrix P* t from (|3"T1) . Lemma (|3.5[) shows that the solution of (1551) is weakly 
closed if the initial data is weakly closed and / = 0. Alternatively, we may use one 
of the implicit Eulerian schemes (l53|) or ((70)) 

(87) (M„ + AtC a -i - AiL a )h* = M Al h*- A * + Aif , 

where L a is the stiffness matrix for either the standard ([53|) or the downwind (fTUj) 
discretization. While these schemes will not preserve the weakly closed property, 
the semi-implicit upwind Eulerian scheme (|72p with algebraic system 

(88) (M„ + AtC CT -i)h* - (M„ + AiL a )h'- At + Ati 
does, due to the interpolation based construction. 

5.2. a-based variational formulation. Here, for simplicity, we assume homoge- 
neous magnetic boundary conditions on dfl. The ansatz e = —Dpa and b = da 
solves Faraday's law ((77)1 . We then multiply Ampere's law with a test form a 1 E T) 1 
and integrate by parts: 

(89) / ft A da' - / (j + /) A a' = 0, Va'eJ^O). 



The material laws ((79")) and (|80l) eliminate j and h and we get the a-based variational 
formulation: Find a <E ^(il) such that: 

(90) / * l± -idaAda' + / * <J D p aAa'= / / A a , Va'eJ" 1 ^). 



The algebraic system for the direct Semi-Lagrangian scheme ([3"2")) then reads as 

(91) (M„ + AtC/i-i y = M ff pi t _ At a*- Ai + Atf , 

while the systems for the implicit and semi-implicit Eulerian schemes are: 

(92) (M CT + AiC M -i + AiL d )a* = M^a'"^ + Atf 
and 

(93) (M CT + AiC M -i)a* = (M ff - AtL d )a*~ A * + Af . 
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Here M CT and C^ 1 are the stiffness matrices for J n * a U} l A % and J Q * fJi -idu)h A 
drjh, uJh,fjh € W 1 (ri/ l ). Ld corresponds either to the standard (j52")l or the upwind 
Eulerian (|69l) discretization of the convection part. 



6. Numerical experiments 

We finally we present a few numerical examples that illustrate the performance 
of the methods derived here. We take fl C R 2 and look at 

dt*uj{t)-£{-l) l d*duj{t) + L p *uj{t) = ip in J] C I 2 , 

(94) i*uj = on ffi, 

w(0) = ujq in fl 

for 1-forms u € .F (fi). In vector proxy notion with u := v. p.(w) this reads 

<9tU + ecurlcurlu + /3(divu) + curl(u x (3) = f in O , 

(95) u • n = on dfl , 

u(0) = u in Q , 

with curl := Rgrad, curl := divR and u x (3 := u T R/3, R = ( !? n ] . We 



v 1 , 

approximate w by discrete 1-forms uih on a triangular mesh. We study the adjoint 
Semi-Lagrangian and Eulerian methods ([55)1 , (1551) , (T7CT|) and (|72|) for this boundary 
value problem. 

Experiment I. In the first experiment we take in problem the domain 
fl = [—1, l] 2 , the divergence free velocity 



(96) (3 = 

and 



(i-x 2 ) 2 (y-y 3 ) 

-(l-y 2 ) 2 {x-x 3 ) 



, n ( sin(7ra;) wsxi-Kii) \ 
cos( 27 rt) ( (1 \ 2 J )(1 A^J 



Then with e = 1 neither the convection nor the diffusion part dominates and we 
expect for all three schemes first order convergence in the mesh size h, if the time 
step At := t — t is of order h. The convergence plots in figures [7] and [5] show the 
L 2 -error at time t = 0.25 and confirm this for different ratios CFL := At ^f^°° . 

h 

Experiment II. We consider the non-divergence free velocity field 



( sin(TTx)(l — y 2 
ysin(Try)(l — x 2 



Again, we obtain first order convergence (see figs [HJ) and fTu) . 

Experiment III. Next, we examine the stability properties of the three different 
schemes for dominating convection, meaning that we choose e in problem (1951) very 
small. While for e = 1 all three schemes are stable for small and larger CFL- 
numbers (see fig. © and ([8])) we encounter this for small e only for the semi- 
Lagrangian and the implicit Eulerian scheme (see fig. (|11|) and (IT2l ). 

As expected in the convection dominated case, one has to use either the semi- 
Lagrangian or fully implicit Eulerian scheme. 

Experiment IV. Another question that we address concerns the stabilization 
for the Eulerian schemes. For the implicit timestepping schemes we need to solve 
in each time step a stationary convection-diffusion problem. It is well-known that 
in the scalar case non-stabilized methods would produce highly oscillating solutions 
for convection-diffusion problems with dominating convection. To study this issue, 
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Adjoint method 



Semi-Lagrangian 
Eulerian Upwind, implicit 
Eulerian Upwind, semi-implicil 
Eulerian Standard, implicit 





10" 
h 



Figure 7. Experiment I: L 2 -error at t = 0.25 for the problem 
with divergence free velocity and small CFL-number CFL = 0.1, 
e = 1. 



Adjoint method 



Semi-Lagrangian 
Eulerian Upwind, implicit 
Eulerian Upwind, semi-implici 
Eulerian Standard, implicit 




10 
h 



Figure 8. Experiment I: L 2 -error at t = 0.25 for the problem 
with divergence free velocity and larger CFL-number CFL = 0.8, 
s = 1. 



we consider the stationary convection diffusion problem for 1-forms and Q C K 

u + e curl cur lu + grad(/3 • u) — R/3curlu = f in Q = [— 1,1] 2 , 
(97) u ■ n = on dfl , 

u(0) = u . 
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Adjoint method 







^M^— Semi -Lag rang ian 

Eulerian Upwind, implicit 

^^^» Eulerian Upwind, semi-implici 
Eulerian Standard, implicit 



















10" 
h 



Figure 9. Experiment II: L 2 -error at t — 0.25 for the problem 
with non-divergence free velocity and small CFL-number CFL = 
0.1, £ = 1. 



Adjoint method 



Semi-Lagrangian 
Eulerian Upwind, implicit 
Eulerian Upwind, semi-implici 
Eulerian Standard, implicit 




10 
h 



Figure 10. Experiment II: L 2 -error at t = 0.25 for the problem 
with non-divergence free velocity and larger CFL-number CFL = 
0.8,e = 1. 



We could either use the standard approximation PO")) or the upwind interpolated 
discrete Lie-derivative ([551 in the discretization of (|9"7)l. Again, we choose 



/ sin(Trx) sin(Try) \ 
U -\(l-a*)(l-y*)) 
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Adjoint method 



Semi-Lagrangian 
Eulerian Upwind, implicit 
Eulerian Upwind, semi-implici : 
Eulerian Standard, implicit 
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FIGURE 11. Experiment III: L 2 -error for different CFL-numbers 
at t = 0.5 for evolution in fixed triangulation with mesh size h m = 
0.11 with e = 10" 3 . 



Adjoint method 



— — Semi-Lagrangian 
-6 — Eulerian Upwind, implicit 
Eulerian Upwind, semi-implici 
Eulerian Standard, implicit 




Figure 12. Experiment III: i 2 -error for different CFL-numbers 
at t = 0.5 for evolution in fixed triangulation with mesh size h m = 
0.11 with e = 10~ 9 . 



and the divergence free velocity (|96|). 

Fig. [13] shows that in the balanced case (e = 1) both the standard and the 
stabilizing upwind scheme yield meaningful solutions. For small e the standard 
method does not seem to converge while the upwind scheme converges (see Fig. [Til) . 
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Figure 13. Experiment IV: i/(curl)-error for stationary problem 
with e = 1. 




Figure 14. Experiment IV: i/(curl)-error for stationary problem 
with e — 10~ 5 . 
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